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Abstract 


Variational models are one of the most efficient techniques for image de- 
noising problems. A variational method refers to the technique of optimizing 
a functional in order to restore appropriate solutions from observed data that 
best fit the original image. This paper proposes to revisit the discrete total 
generalized variation (TGV) image denoising problem by redefining the op- 
erations via the inclusion of a diagonal term to reduce the staircasing effect, 
which is the patchy artifacts usually observed in slanted regions of the image. 
We propose to add an oblique scheme in discretization operators, which we 
claim is aware of the alleviation of the staircasing effect superior to the con- 
ventional TGV method. Numerical experiments are carried out by using the 
primal-dual algorithm, and numerous real-world examples are conducted to 
confirm that the new proposed method achieves higher quality in terms of rel- 
ative error and the peak signal to noise ratio compared with the conventional 
TGV method. 
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1 Introduction 


Digital image processing (DIP) deals with performing operations on digital 
images. A digital image is a numerical representation of a physical scene, 
which is composed of a finite number of pixels. Digital images are produced 
by means of imaging machines that cover the electromagnetic spectrum. Syn- 
thetic images, electron microscopy images, and ultra-sound images are exam- 
ples of digital images. Digital images were first used in the newspaper indus- 
try in the 1920s. These digital images were produced from a coded tape by 
a telegraph printer. The field of DIP is enriched with various applications, 
including image restoration [3, 17,23, 25], artistic effects [19], medical visual- 
ization [9, 27,28, 30], industrial inspection [22], law enforcement [32,33], and 
so on. Image restoration is one of the most widespread applications of DIP 
techniques that implements processes on digital images in order to estimate 
the original image from the corrupted one. The image distortion is caused 
due to different types of noise, such as Gaussian noise, white noise, salt and 
pepper noise, and speckle noise. In recent decades, variational approaches 
have been used as an efficient tool for image denoising problems. 


A variational model is an optimization problem in which the criterion 
is defined as a functional (energy), which consists of a regularization term 
and a data fidelity term. Total variation (TV) regularization is a variational 
model that uses total variation as a regularization term. TV regularization 
was first proposed by Rudin, Osher, and Fatemi (ROF model) for imaging 
problems; see [26]. The ROF method is edge-preserving and has a fast nu- 
merical algorithm. Many different papers have shown the efficiency of TV 
minimization for image restoration [1,4,5,8, 10,11, 13-15, 18,20, 21,24, 29,31]. 
The TV regularization has been widely used in various applications, such as 
image deblurring, inpainting, image zooming, segmentation problems, inter- 
polation, spectral extrapolation, and stereovision. In these methods, the TV 
semi-nom is defined as 


TV (u) = sup {| u div v dx | vEC(A,R"), |lvIloo < i, 
Q 


where u is a function defined on a bounded region Q C R”. TV based 
regularization models have been proved to be efficient in image denoising 
problems. However, these models suffer from the staircasing effect, which 
appears as undesired patchy artifacts in slanted regions (see Figure 1). The 
total generalized variation (TGV) regularization model [6] is one technique 
to overcome this shortcoming, which acts as a regularization functional that 
incorporates higher-order derivatives and regularizes independently on vari- 
ous regularity levels. The main idea of TGV is a generalization of TV, which 
is defined as 
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Figure 1: The implementation results of 1500 iterations of the ROF TV model [26]. From 
left to right: The reference image, the noisy image, and the restored image using the ROF 
TV. The original image has been corrupted by additive white Gaussian noise of standard 
deviation 0.18 and the regularization parameter is set to A = 0.17. The ROF denoising 


model leads to staircasing effect, which is observed as patchy artifacts. 


TGV¥ (u) =sup {| 


a 
1=0,...,k-1}, 


u div’v dx | vy € CEQ, Sym*(R")), ||div'v|loo < AL, 


where k > 1, n > 1, and %; are fixed positive parameters, and Sym*(R”) 
denotes the space of symmetric tensors of order k with arguments in R”. 
This paper proposes to revisit the TGV image denoising model by redefin- 
ing the gradient operations via the inclusion of diagonal terms to reduce the 
staircasing effect. We propose to add an oblique scheme in classical image 
derivatives discretization, which we claim is aware of the alleviation of the 
staircasing effect superior to the conventional TGV method. Numerical ex- 
periments are carried out using the primal-dual algorithm, and numerous 
real-world experiments are conducted to confirm the effectiveness of the new 
approach. 


The remainder of this paper is organized as follows. First, a brief ex- 
planation of some essential concepts regarding the TV and T’'GV schemes is 
presented in section 2. In section 3, the new proposed regularization model 
for image denoising problems and the corresponding theoretical results are 
presented. Section 4 is devoted to numerical experiments and comparisons 
that demonstrate the efficiency of the proposed method. Finally, Section 5 
contains some concluding remarks. 


2 Background 


In this section, we present a brief review on essential concepts of TV and 
TGV approaches. 
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2.1 TV concept 


Assume that u is a function defined on a bounded region 2 C R”. The 
function wu is said of bounded variation (BV function) if it is integrable and 
there exists a Radon measure Du such that 


[uvaiwvae=— fv Dude, 
Q Q 


for ally € C}(R”,R”), where C}(R”, R”) is the space of continuously differ- 
entiable vector functions v of compact support and Du is the distributional 
(weak) derivative of u. The TV seminorm of u is defined as 


TV(u) := sup | u div v dx | v € 02(Q,R"), |\vIloo < i} = i, |Duldz, 
Q 


where |.| is the Euclidean norm. In the case that u is a smooth function, we 
have Du = Vu; therefore TV (wu) is the integral of its gradient magnitude 


TV(u) = f [ula 


2.2 TGV concept 


In this subsection, we review the essential concepts of TGV regularization 
from [6]. 

Definition 1. Let k > 1, Qc R” be a bounded region, and Xo,..., Ax—1 be 
fixed positive parameters. For u € Lj,,.(), the total generalized variation of 


order k with weight A € R* is defined by the functional 


TGV (u) = sup { | 


u div'v dx | vy € CF(2, Sym*(R")), ||div'V|loo < 
Q 


SG ngk tts 


where \ = (Ao,---; Ax—1), C*(Q, Sym*(R"”)) is the space of k-times contin- 
uously differentiable functions with compact support from Q to Sym*(R”) 
and Sym*(R”) is the vector space of symmetric k-tensors defined as 
Sym*(R”) = {n >R°x-- xR" OR | 7 is k-linear and symmetric}. 
emC—~__—__ 


k times 


For the case when k = 2, the TGV¥ functional is a special case of TGV," 
functional, which is defined as 


TGV; (u) = sup { | u div?v dx | vy € C2(Q, Sym?(R”)), ||V|loo < Ao, 
Q 
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IIdiv v|loo < Ax}, (1) 
where \ = (Ao, Az) and 
OV 4 ay O14; 
d i= aj d = Ut a7 
(ee 2 an? Le Ou? 1 BuiOx 


Remark 1. The Euclidean space R™*4?2 ig equipped with the inner product 


N, No 
N,{XN: 
(U, 8) RN xNo = ) ) UijSi,j, USE Roe. 


i=1 j=l 


For u € R™!*%? | the second order TGV semi-norm (1) can be discretized 
as 


TGV (u) = max { (u, div?) prise ge ye, p= & i ; 


V2 V22 


vig E RN? (i,j = 1,2), |[Vlloo < Ao, |Idiv vlloo Sa}, (2) 


where \ = (Ag, Ai). The discrete version of the infinity norm for the vector 
field z, where z = div v, is defined as 


1/2 
la = sn | ( Sac? | real 


Moreover, the discrete version of the infinity norm for matrix v is defined as 


n 1/2 
ila = sue {( Sve? +2 ve? } | reat 


i<j 


3 The new proposed method 


In this section, we present a new variant of the TGV model for the alleviation 
of the staircasing effect in image denoising problems by redefining the oper- 
ators via the inclusion of a diagonal term. In the conventional TGV model 
presented in [6], the discretization operators are based on finite differences 
in the direction of the horizontal and vertical axes. Here, we reconstruct the 
discretization operators in [6] via the inclusion of a diagonal term. For this 
purpose, we introduce the versions of the discretization operators in parts (a- 
e). The indexes x, y, and o indicate that the corresponding finite-differences 
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are established in the direction of the horizontal, vertical, and diagonal axis, 
respectively. 
Ni x No : . & : Ni x No 
a. Forue R , we define the discretization operator 6: R — 
(Reyer as 


LS Wi41,7 7 Ui, j if 1<i<M,, 
i> hee eee 


+.)  _ J Uijt1 — Uij if 1 <j < No, 
(His = { 0 if j= Np, 


(Otu)ij = UWit1,j+l — Wij if L<i<N,,1<j< No, 
0 Wi.3 0 if i= M,j = No. 


Pi 
b. For p= | po |, p; € R™**? (i = 1, 2,3), the discretization operator 


P3 
¢: (R3)%1*N2 —, R%1*N2 is defined as 


C(p) = 9; (p1) + Oy (p2) + O5 (ps), 


where 
@ieg-— Onaga lt tL 4e N4, 
(O5p1)4 = (pi)sg if t=1, 
—(piji-1,y if += M1, 
(p2)i,j — (Pa)ig-1 if 1<7< Na, 
(0; p2)ig = (p2)i,g if 7 =1, 
—(po)ig-1 if 7 = No, 
(p3)i,j — (Pa)i-1,j-1 if 1<ti<M,1<j< No, 
(OF p3)ig = 
eg LiSly=1. 
Pl 
c. For p= | po} , p; €e R™*%? (i = 1,2,3), the discretization operator 
P3 


&: (R°)M1*N2 —, (R°)1*42 is defined as 


A new regularization term based on second order total generalized variation ... 147 


d; (p1) 5(O; (pa) + Oy (p1)) 5 (Ox (ps) + 95 (p1)) 
E(p) = | 3(Oz (p2) + OF (p1)) d; (p2) 5(O, (ps) + 5 (p2)) | 
3(Oz (ps) + 85 (p1)) 3 (Oy (ps) + OF (p2)) a; (ps) 


where 0; , 0, , and 0; are defined as in part (b). 


41 12 13 
d. For v = | 42 v22 ¥23] , Vag € RN1*N2 (i,7 = 1,2,3), the discretiza- 
113 V23 33 
tion operator ¢ : (R°)N1*N2 —» (R®)%1*2 is defined as 


Oz (vir) + OF (Viz) + OF (13) 
C(v) = | OF (2) + OF (V22) + OG (Vea) | ; 
Oz (vis) + Oy (v23) + O5 (V33) 


where O07, 0, and OF are defined as in part (a). 


141 12 V13 
N,x Noe eh o> 8 . . 
e. For vy = | 42 V22 423] , Wy CE R (i,7 = 1,2,3), the discretiza- 
113 V23 V33 
tion operator ¢? : (R°)%1*N2 —; RN'*N2 is defined as 


C7(v) = 0; OF (v1) + OF OF (v22) + OF OF (v33) + OF OF (viz) + OF OF (13) 


+0; OF (12) + O, OF (ves) + 0; OF (113) + O, OF (v23). 
where Of , O°, OF, 0; , O,, and 05 are defined as in parts (a) and (b). 


The next step, we aim to solve the following variational image denoising 
problem 
min $||u— uo||3 + La(u), (3) 
u 


where up € R%'!*? is the noisy image, u € R™:*? is the image to be 
reconstructed, X = (Ag, A1) (Ao, A1 are the regularization parameters), and 
E(u) is the new proposed regularization functional, which is defined as 


11 12 13 


L)(u) = max {(u , Cv) pmixny, |v € (R°)M*N2, vy = | 412 Vo2 V23 | , 
113 V23 V33 


Vig E RNXN (i,j = 1,2,3), [[Vlloo S ro, II¢(v)Iloo < ai}. (4) 
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3.1 Theoretical results 


We apply the over-relaxed Chambolle—Pock algorithm described in [16] for 
solving (3). Since this algorithm solves jointly the primal and dual formu- 
lations of minimization problem (3), we need to obtain the Fenchel dual 
formulation of functional (4). The analytical process for establishing this 
formulation of (4) is stated in Theorem 2, in which we follow the steps of [7] 
for its proof. Before studying this theorem, the reader needs to be familiar 
with the concept of Legendre—Fenchel duality and Fenchel duality theorem 
from [2]. 


Definition 2. Given some convex, proper, and lower semi-continuous func- 
tion f(p) defined for p € H, where H is a Hilbert space with inner product 
(.,+) 7» its Legendre-Fenchel dual function is defined as 


f*(q) = max {(p, 9) 7 — f(p) | p © A} 
for all q € H. 


Theorem 1. Assume that X and Y are real Banach spaces, that f, : X —> 
(—oo, +00) and fo : Y —> (—oo, +00) are proper, conver, and lower semi- 
continuous functions, and that A: X —>Y is a linear continuous operator. 
If there exists 19 € X such that fy(ap) < co and f2 is continuous at Axo, 
then 


min {fi(w) + fo(A(a)) | « € X} = max {—f(y") — fi(-A*y") | ye YF. 


Pi Wi 
Remark 2. For p= | po | € (R®)%2*4? and w = | we | € (R°)™*%2, the 
P3 W3 


Euclidean space (R°®)%1*%? is equipped with the inner product 


Ni No 


(D,W) (R38) Ni x2 = S> S (ida (wr)ij + (p2)i,j(wa)i,g + (ps)i,j (ws )a7- 


i=1 j=1 


Theorem 2. The discrete second order total generalized variation functional 
(4) is equivalent to the following Fenchel dual formulation 


Ly(u) = min Apol|E(p)||1 + A1||6(u) — pla, 
pe (RE)NixN2 | 


where € and 6 are defined as in parts (c) and (a), respectively. 


Proof. We prove in the lines of [7]. Based on (4), we have 
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Ly(u) = max {(u, C20) paix |v ERMAN, [Iya < do; [[G(Plloo < Ar} 


= max { (u,62(Y) pews xna — II, [1oo<A0}) — Lelhtloosarp(6(0)) | » € (ROY *N2 f 


ve (R°)N1*Na i. 


= = min { 1,\) J10<ro}) + TE) ).|Joo< ar} (S()) _ (u, Cv) pu, x No 
where ea 
_ 0 if |lv oo S Ao; 
Li) \Jooso}(¥) = { if |IVlloo > Ao, 


Hasan) = {9 it Hewole > 


We choose 


fil¥) =I flecsro¥)s f2(C(Y)) = Le [osa3 (CY) — (us P()) para xine - 


Based on the principles of Theorem 1, it follows that 


Ly(u) = —min { filv) + fo(¢(v)) | vy € (R°)NixN2 } 


= min { fi(-€(P)) + £8@)| p  (RAYM*™ J, 
where p = C(v). Based on the duality principle of definition 2, it yields 
fi (—€(p)) = Aoll€(P)||1, and 
f3(p) = max {(p, w) aqsyvixve — falw) | we (RAYN*N2} 
= max {(p, W) (Ra) Ni xo _ Teh) joo<ar} (WY) 1 (u, ¢(w)) Rm xNo | WE (R*) 
Since ¢* = —d, it yields 
f3(p) = max {(4(u), w) aasyai xno ~~ (P, W) (Rs) M1 x Ne aes CA, WE (RAPOSs} 
Ilelleo < Aa, w © (R)N1*N2}, 


Nene. 


= max {(4(u) — p, w)(Rsywi x2 
Choosing k = 6(u) — p, it follows that 


2 No M 


f3(p) = max { >> Yweli, 3) keli,3) | twlloo < Ar, w € (RP) *™} 


s=1 j=1 i=l 


Ai ||d(u) — p|fi- 


3.2 Primal-dual algorithm 


This section contains the primal-dual algorithm described in [16] for solving 
(3). Primal-dual methods apply proximity operators, which can be defined for 
proper, lower semi-continuous, convex, and extended real-valued functions. 
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The followings are the proximity operators, which are used in this algorithm, 


U+ TUO _ Pp 
prot, rp, (u) = adage Pees (p) =p- [pl ’ 
max(——, 1) 
TAY 
v 
PTOLgg( j= Iv] 
ax(—,1 
m, Gr ) 


Algorithm 1 to solve (3) 


1. Set &k = 0, choose parameters 7, 0, p, and the initial estimates 
Pn) E Beene. p) c (Reyne, p() € (R°)MixNe 

2. Calculate ut), p+) and v(t) using the following equations: 

Py := prozrr,(u® — ¢(r¢(v))), 

Py := prowrr,(p® + 7¢6(v)), 

Ps = protag(v'") +0 €(5(2 pr — u) ~ (2 pp — p™))), 

ult) c= ul) + p (py) — ul), 

pe =o" +p ae), 

RTL) -— yk) 4 9 (p3 — yi*)), 

3. Stop or set k :=k+1 and go back to step 2. 


4 Numerical experiments 


In this section, we test the performance of the proposed method on several 
sample images to remove noise (see Figure 2). Each image is corrupted by ad- 
ditive white Gaussian noise of standard deviation 0.18. We compare the new 
proposed method with the anisotropic TV, the isotropic TV, the upwind TV 
n [12], and the TGV method in [6] (see Figures 3,4,5,6,7,8,9,10,11,12,13,14 
for the restored images). For the selection of the optimal regularization pa- 
rameter, the algorithms corresponding to the anisotropic TV, the isotropic 
TV, the upwind TV, the TGV method, and the proposed method are imple- 
mented with many different choices for A, and the \ value corresponding to 
the best peak signal to noise ratio (PSNR) (least relative error) is chosen as 
the optimal X value. 

Using the anisotropic TV, it performs well in removing noise but small 
details become unclear, and this TV scheme suffers from the staircasing ef- 
fect, which appears as undesired patchy artifacts in slanted regions. The 
upwind TV has nice performance in preserving details but it has a remark- 
able drawback because some small white particles of noise will be remained 
in restored images, which indicates that the upwind TV is not capable of 
removing white noise. Moreover, the upwind TV suffers from the staircasing 
effect. The isotropic TV performs to some extent better than the anisotropic 


A new regularization term based on second order total generalized variation ... 151 


TV and the upwind TV in noise removal but still some details are not recon- 
structed during the denoising process, and this TV scheme is not capable of 
alleviating the staircasing effect. Using the TGV method, it still suffers from 
the remaining of the staircasing effect, which indicates that it is not powerful 
enough to handle the severe noise. The proposed method outperforms the 
anisotropic TV, the isotropic TV, the upwind TV, and the TGV method both 
in the reconstruction of fine structures and the elimination of the staircasing 
effect. The numerical experiments illustrate that the new proposed method 
achieves higher quality in terms of PSNR and relative error (see Table 1). 
For example in the case of bird image, Table 1 illustrates that the proposed 
method achieves a PSNR value, which is about 0.187585 higher than the 
PSNR value of TGV model. If we notice the other quantities in Table 1, 
we observe that the PSNR value corresponding to the TGV model is about 
0.129112 higher than the PSNR value corresponding to isotropic TV, and this 
PSNR. difference is less than the PSNR difference of the proposed method 
and the TGV method. 

The iteration numbers of optimization algorithm for both bird and bike 
images in each model is 1500. Figures 15,16,17, and 18 illustrate the PSNR 
and the relative error of various methods versus iteration numbers. For ex- 
ample, in the case of bike image, Figures 17 and 18 illustrate that after 600 
iterations, the proposed method has smaller PSNR and higher relative er- 
ror in comparison with the TGV model (TGV: PSNR=23.038355, Relative 
error=0.092436, proposed: PSNR=22.980077, Relative error=0.093058). If 
we increase the number of iterations to 1500, the TGV model achieves better 
quantities (TGV: PSNR=23.060928, Relative error=0.092196) in compari- 
son with 600 iterations, and the proposed model achieves the best PSNR 
value (proposed: PSN R = 23.080060, Relative error=0.091993) in compari- 
son with other methods. Even if we increase the iteration numbers again, very 
few changes are observed, and the proposed method will still have the best 
PSNR value and the least relative error in comparison with other methods. 
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Figure 2: The test images used in our numerical experiments. Left: bird (490674 in 
experiments); right: bike (640x512 in experiments). 
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Figure 5: From left to right: restored bird images by TGV and the new proposed method. 
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Figure 6: Zoomed-in regions of bird image. From left to right: noisy image, restored by 
anisotropic TV. 


Figure 7: Zoomed-in regions of bird image. From left to right: restored by isotropic TV 
and upwind TV. 


Figure 8: Zoomed-in regions of bird image. From left to right: restored by TGV and the 
new proposed method. 
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Figure 9: From left to right: noisy bike image, restored bike image by anisotropic TV. 


Figure 10: From left to right: restored bike images by isotropic TV and upwind TV. 
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Figure 11: From left to right: restored bike images by TGV and the new proposed 
method. 


Figure 12: Zoomed-in regions of bike image. From left to right: noisy image, restored 
by anisotropic TV. 
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Figure 13: Zoomed-in regions of bike image. From left to right: restored by isotropic 
TV and upwind TV. 


al 


Figure 14: Zoomed-in regions of bike image. From left to right: restored by TGV and 
the new proposed method. 
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Figure 15: The sequence of relative error of bird image versus iteration numbers. 
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Figure 16: The sequence of PSNR of bird image versus iteration numbers. 
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Figure 17: The sequence of relative error of bike image versus iteration numbers. 
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Figure 18: The sequence of PSNR of bike image versus iteration numbers. 
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Table 1: Performance comparison of restoration results in terms of relative error and 
the PSNR. The new proposed method achieves higher quality in terms of the PSNR and 


relative error. 


Image Method Iterations Optimal A Relative error PSNR 
Bird Anisotropic TV 1500 0.14 0.064262 27.313069 
Bird Isotropic TV 1500 0.17 0.063141 27.466009 
Bird Upwind TV 1500 0.21 0.065420 27.158004 
Bird TGV 1500 (0.3, 0.16) 0.062209 27.595121 
Bird proposed 1500 (0.3, 0.13) 0.060880 27.782706 
Bike Anisotropic TV 1500 0.11 0.092527 23.029824 
Bike Isotropic TV 1500 0.13 0.092511 23.031338 
Bike Upwind TV 1500 0.155 0.097285 22.594246 
Bike TGV 1500 (0.21, 0.13) 0.092196 23.060928 
Bike proposed 1500 (0.14, 0.1) 0.091993 23.080060 


5 Conclusion 


This paper proposes to revisit the discreteTGV image denoising problem by 
redefining the operations via the inclusion of a diagonal term to reduce the 
staircasing effect, which is the patchy artifacts usually observed in slanted 
regions of the image. Numerical experiments confirm that the new proposed 
method achieves higher quality in terms of relative error and the PSNR com- 
pared with the conventional TGV method. The more direction analysis of 
first-order-derivative by using more oblique terms is considered as our future 
researches. 
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